rm(list = ls())
options(width=100, device = 'quartz')

###################
## Assemble the Data 
####################
library(foreign)
library(plyr); library(dplyr)
library(coin)
library(rdrobust)
library(rms)
library(rdd)
library(ggplot2)
library(reshape2)
library(xtable)
library(MASS)
library(plm) ## for panel models
library(stargazer) ## for nice latex tables
library(lmtest) ## for coeftest
library(clusterSEs) ## for bootstrapped SEs

myRecode <- function(var, recode.ls, as.factor.result=TRUE, ...) {
  options(useFancyQuotes = FALSE)
  recode.ls <- sapply(recode.ls, function(x) paste(sQuote(x), collapse='='))
  recodes <- paste(recode.ls, collapse=';')
  print(recodes)
  revar <- Recode(var, recodes, as.factor.result=TRUE, ...)
  options(useFancyQuotes = TRUE)
  revar
}
rols <- function (..., data, subset=TRUE, cluster=1:nrow(data),
                  method="huber") {
  mod <- ols(..., data=data[subset, ], x=TRUE)
  robcov(mod, cluster=cluster[subset], method=method)
}
PasteEval <- function(..., print = FALSE) {
  txt <- paste(..., sep = "")
  if (print) { cat("\n", txt, "\n") }
  eval(parse(text = txt), envir = parent.frame())
}
FileName <- function (path=NULL, name=NULL, ext=NULL, replace=FALSE) {
  p <- paste(path, collapse = "")
  n <- paste(name, collapse = "")
  e <- paste(".", ext, sep = "")
  d <- format(Sys.Date(), "%y%m%d") 
  for (l in seq_along(letters)) {
    if (l > 1) old_file_name <- file_name
    file_name <- paste0(p, d, paste(n, letters[l], sep="-"), e)
    if (!any(grepl(file_name, list.files()))) {
      if (replace) file_name <- ifelse(l > 1, old_file_name, file_name)
      break
    }
  }
  cat("\nFile name:", file_name, "\n\n")
  return(file_name)
}

mysw <- function (...) {
  setwd(paste(..., sep = "/"))
}
mypdf <- function(name, ...) {
  date <- format(Sys.Date(), "%y%m%d")
  pdf(paste(paste(name, collapse=""), date, ".pdf", sep=""), ...)
}
tri <- function (x, h, c=0) pmax(0, 1 - abs((x - c) / h))
POtoSouth11 <- function (pos, dta=st.info) {
  dta$South11[match(pos, dta$POAbrv)]
}



mysw(home.dir)

st.info <- read.csv("StateCodes.csv")


load(file= "party_control_data_161121.RData")

data.use <- data %>%
    plyr::arrange(abb, year) %>%
        dplyr::filter(!is.na(abb) & year >= 1936)

table(data.use$year)

################
#### Setup RDD Variables
################

data.use <- data.use %>%
    group_by(
        abb
    ) %>% mutate(
        South11 = POtoSouth11(abb),
        DemGov = gov_party,
        DemPropGov0 = demprop2 - 1/2,
        DemMarginGov = DemPropGov0 * 2,
        DemWinGov = as.integer(DemPropGov0 > 0),
        DemSeatShareHouse = hs_dem_per_2pty / 100,
        DemSeatShareSenate = sen_dem_per_2pty / 100,
        DemControlHouse = hs_dem_control,
        DemControlSenate = sen_dem_control,
        GovLib = Policy,
        GovLibL1 = lag(GovLib, 1),
        GovLibL2 = lag(GovLib, 2),
        GovLibL3 = lag(GovLib, 3),
        GovLibL4 = lag(GovLib, 4),        
        GovLibL12 = (GovLibL1 + GovLibL2) / 2,
        GovLibP1 = lead(GovLib, 1),
        GovLibP2 = lead(GovLib, 2),
        GovLibP3 = lead(GovLib, 3),
        GovLibP4 = lead(GovLib, 4),
        GovLibP12 = (GovLibP1 + GovLibP2) / 2,
        GovLibD1 = GovLibP1 - GovLib, ## delta 1
        GovLibD2 = GovLibP2 - GovLib,  ## delta 2
        GovLibD3 = GovLibP3 - GovLib,  ## delta 3
        GovLibD4 = GovLibP4 - GovLib,  ## delta 4
        GovLibD12 = (GovLibD1 + GovLibD2)/2, ## ave of deltas 1-2
        GovLibDL1 = GovLib - GovLibL1 ## change, t relative to t-1
    ) %>% ungroup()


res.fe <- resid(lm(GovLib ~ abb + factor(year), data=data.use))
res.l1 <- resid(lm(GovLib ~ Policy_L1, data=data.use))
res.l2 <- resid(lm(GovLib ~ Policy_L2, data=data.use))
res.fe.l1 <- resid(lm(GovLib ~ abb + factor(year) + GovLibL1, data=data.use))
res.fe.l2 <- resid(lm(GovLib ~ abb + factor(year) + GovLibL2, data=data.use))

data.use$GovLib.fe0 <- NA
data.use$GovLib.fe0[as.integer(names(res.fe))] <- res.fe
data.use$GovLib.fe1 <- NA
data.use$GovLib.fe1[as.integer(names(res.fe.l1))] <- res.fe.l1
data.use$GovLib.fe2 <- NA
data.use$GovLib.fe2[as.integer(names(res.fe.l2))] <- res.fe.l2
data.use$GovLib.1 <- NA
data.use$GovLib.1[as.integer(names(res.l1))] <- res.l1
data.use$GovLib.2 <- NA
data.use$GovLib.2[as.integer(names(res.l2))] <- res.l2

data.use <- group_by(data.use, abb) %>%
    mutate(GovLib.1L1 = lag(GovLib.1, 1), ## lag of first-difference
           GovLib.1L2 = lag(GovLib.1, 2),
           GovLib.2L1 = lag(GovLib.2, 1), 
           GovLib.2L2 = lag(GovLib.2, 2),
           GovLib.fe0L1 = lag(GovLib.fe0, 1),
           GovLib.fe0L2 = lag(GovLib.fe0, 2),
           GovLib.fe1L1 = lag(GovLib.fe1, 1),
           GovLib.fe1L2 = lag(GovLib.fe1, 2),
           GovLib.fe2L1 = lag(GovLib.fe2, 1),
           GovLib.fe2L2 = lag(GovLib.fe2, 2),
           GovLib.1P1 = lead(GovLib.1, 1), ## lead of first-difference
           GovLib.1P2 = lead(GovLib.1, 2),
           GovLib.2P1 = lead(GovLib.2, 1), 
           GovLib.2P2 = lead(GovLib.2, 2),
           GovLib.fe0P1 = lead(GovLib.fe0, 1),
           GovLib.fe0P2 = lead(GovLib.fe0, 2),
           GovLib.fe1P1 = lead(GovLib.fe1, 1),
           GovLib.fe1P2 = lead(GovLib.fe1, 2),
           GovLib.fe2P1 = lead(GovLib.fe2, 1),
           GovLib.fe2P2 = lead(GovLib.fe2, 2)
           )

tail(dplyr::select(data.use, c(year, GovLib, GovLib.fe0, GovLib.fe0L1, GovLib.fe0L2,
                              GovLib.fe0P1, GovLib.fe0P2)))



################################################################################
#### ESTIMATES #################################################################
################################################################################

####
### PLACEBO OUTCOMES
#####

### PARTY CONTROL
with(data.use, rdrobust(y=DemGov, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=DemSeatShareHouse, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=DemControlHouse, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=DemSeatShareSenate, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=DemControlSenate, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=DemControl, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

### POLICY
## Contemporaneous
with(data.use, rdrobust(y=GovLib, x=DemMarginGov,
                    bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib - GovLibL1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib - GovLibL2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib.1, x=DemMarginGov,
                        bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib.fe0, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> **borderline**
with(data.use, rdrobust(y=GovLib.fe1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib.fe2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

### Lag 1
with(data.use, rdrobust(y=GovLibL1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLibL1 - GovLibL2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib.1L1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib.2L1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib.fe0L1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> **borderline**
with(data.use, rdrobust(y=GovLib.fe1L1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data.use, rdrobust(y=GovLib.fe2L1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

placebo.df <- subset(data.use,, c(DemGov, DemControlHouse, DemSeatShareHouse,
                                  DemControlSenate, DemSeatShareSenate,
                                  GovLib, GovLibDL1))

placebo.res <- llply(placebo.df, function (placebo) {
  rdrobust(y = placebo, x = data.use$DemMarginGov,
           bwselect="mserd", all=TRUE, kernel='tri')
})

placebo.tab <- ldply(placebo.res, function (res) {
  data.frame(BW = res$bws[1, "left"],
             Est = res$coef["Conventional", 1],
             LoCI = res$ci["Robust", "CI Lower"],
             HiCI = res$ci["Robust", "CI Upper"],
             PV = res$pv["Robust", "P>|z|"])
})

rownames(placebo.tab) <- placebo.tab[, 1]
placebo.tab <- placebo.tab[, -1]
xtable(placebo.tab)



##################
### ESTIMATES
##################

#### Main estimates for paper:
## Change t + 1
with(data.use, rdrobust(y = GovLibD1, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
## Change t + 2
with(data.use, rdrobust(y = GovLibD2, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))

#### Other estimates:
##> +
## Lead 1
with(data.use, rdrobust(y=GovLibP1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data.use, rdrobust(y=GovLib.1P1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data.use, rdrobust(y=GovLib.fe0P1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + **
with(data.use, rdrobust(y=GovLib.fe1P1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + (similar to FE)
with(data.use, rdrobust(y=GovLib.fe2P1, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))

##> + *
## Lead 2
with(data.use, rdrobust(y=GovLib.fe0P2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + ***
with(data.use, rdrobust(y=GovLib.2P2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + *
with(data.use, rdrobust(y=GovLib.fe1P2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + *
with(data.use, rdrobust(y=GovLib.fe2P2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + ** (similar to FE)

## Average of Lead 1 and 2
with(data.use, rdrobust(y=GovLibP12, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data.use, rdrobust(y=(GovLib.1P1 + GovLib.2P2) / 2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data.use, rdrobust(y=(GovLib.fe0P1 + GovLib.fe0P2) / 2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + **
with(data.use, rdrobust(y=(GovLib.fe1P1 + GovLib.fe1P2) / 2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + **
with(data.use, rdrobust(y=(GovLib.fe2P1 + GovLib.fe2P2) / 2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data.use, rdrobust(y=(GovLib.fe1P1 + GovLib.fe2P2) / 2, x=DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + *
## average of change t + 1 and t + 2
with(data.use, rdrobust(y = GovLibD12, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + 

## PRE-1980
with(subset(data.use, year <= 1980),
            rdrobust(y = GovLibD1, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
with(subset(data.use, year <= 1980),
            rdrobust(y = GovLib.fe1P1, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
with(subset(data.use, year <= 1980),
            rdrobust(y = GovLib.fe2P2, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
## POST-1980
with(subset(data.use, year > 1980),
            rdrobust(y = GovLibD1, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
with(subset(data.use, year > 1980),
            rdrobust(y = GovLib.fe1P1, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))
with(subset(data.use, year > 1980),
            rdrobust(y = GovLib.fe2P2, x = DemMarginGov,
                     bwselect="mserd", all=TRUE, kernel='tri'))

################################################################################
#### PLOTS #####################################################################
################################################################################

getRD <- function (rdr.out, which.get = c('Conventional', 'Robust')) {
  rdr.out$tabl3.str[which.get, ]
}
myRD <- function (ynm, keep = TRUE, data = data.use, xnm = "DemMarginGov", c = 0,
                  get = TRUE, level = 95, ...) {
  rdr <- rdrobust(y = data.use[[ynm]], x = data.use[[xnm]], subset = keep,
                  c = c, bwselect = "mserd", all = TRUE, kernel = 'tri',
                  level = level)
  if (get) {
    getRD(rdr, ...)
  } else {
    rdr
  }
}

n.leads <- 4
(yr.brks <- c(min(data.use$year), 1967, 1989, max(data.use$year)))
data.use <- mutate(data.use,
                   Era = cut(year, yr.brks, include.lowest = TRUE, dig.lab = 4))
yr.subs <- llply(levels(data.use$Era), function (e) {
  sort(unique(data.use$year[data.use$Era == e]))
})
names(yr.subs) <- laply(yr.subs, function (yrs) paste0(min(yrs), "-", max(yrs)))
yr.subs <- c("All Years" = list(sort(unique(data.use$year[!is.na(data.use$Era)]))),
             yr.subs)

(dvs <- paste0("GovLibD", 1:n.leads))
rd.ls <- llply(yr.subs, function (sub) {
  llply(dvs, function (dv) {
    myRD(dv, keep = data.use$year %in% sub, level=95)
  })
})

rd.cast <- dcast(melt(rd.ls), L1 + L2 ~ Var2 + Var1)
rd.cast[, -c(1:2)] <- sapply(rd.cast[, -c(1:2)], as.numeric)
rd.cast <- rd.cast %>% mutate(
    Lead = factor(L2, labels = 1:n.leads),
    L1 = factor(L1),
    Era = factor(L1,
        levels=c(levels(L1)[nlevels(L1)], levels(L1)[-c(nlevels(L1))]))
)
levels(rd.cast$Era)[nlevels(rd.cast$Era)] <-
  gsub(max(data.use$year), max(data.use$year) - 1,
       levels(rd.cast$Era)[nlevels(rd.cast$Era)])
summary(rd.cast)


(yr.brks <- c(min(data.use$year), 1989, max(data.use$year)))
data.use <- mutate(data.use,
                   Era = cut(year, yr.brks, include.lowest = TRUE, dig.lab = 4))
yr.subs <- llply(levels(data.use$Era), function (e) {
  sort(unique(data.use$year[data.use$Era == e]))
})
names(yr.subs) <- laply(yr.subs, function (yrs) paste0(min(yrs), "-", max(yrs)))
yr.subs <- c("All Years" = list(sort(unique(data.use$year[!is.na(data.use$Era)]))),
             yr.subs)

(dvs <- paste0("GovLibD", 1:n.leads))
rd.ls <- llply(yr.subs, function (sub) {
  llply(dvs, function (dv) {
    myRD(dv, keep = data.use$year %in% sub, level=95)
  })
})

rd.cast <- dcast(melt(rd.ls), L1 + L2 ~ Var2 + Var1)
rd.cast[, -c(1:2)] <- sapply(rd.cast[, -c(1:2)], as.numeric)
rd.cast <- rd.cast %>% mutate(
    Lead = factor(L2, labels = 1:n.leads),
    L1 = factor(L1),
    Era = factor(L1,
        levels=c(levels(L1)[nlevels(L1)], levels(L1)[-c(nlevels(L1))]))
)
levels(rd.cast$Era)[nlevels(rd.cast$Era)] <-
  gsub(max(data.use$year), max(data.use$year) - 1,
       levels(rd.cast$Era)[nlevels(rd.cast$Era)])
summary(rd.cast)




mypdf("Figure3_GovRDAll-", width=5, height=3)
(ggplot(subset(rd.cast, Era == "All Years"))
 ## + facet_wrap(~Era, ncol = 4)
 + geom_pointrange(aes(x = Lead, y = Coef_Conventional, ymin = `CI Lower_Robust`,
                       ymax = `CI Upper_Robust`))
 + geom_hline(yintercept = 0, linetype = 'dotted')
 + ylab('RD Effect on Policy Change\n(Robust 95% CI)')
 + xlab("Years after Election")
 + ggtitle("Governor")
 + theme_bw()
 )
dev.off() 


### WINDOWS
yr.window <- 20
window.mids <- seq.int(min(data.use$year) + yr.window / 2,
                       max(data.use$year) - yr.window / 2 - 1)

window.ls <- llply(window.mids, function (mid) {
  seq.int(mid - yr.window / 2, mid + yr.window / 2)
})

rd.win.ls <- llply(window.ls, function (win) {
  myRD("GovLibD1", keep = data.use$year %in% win, level=95)
})
names(rd.win.ls) <- window.mids

rd.win.cast <- melt(rd.win.ls) %>% dcast(L1 ~ Var2 + Var1)
rd.win.cast <- sapply(rd.win.cast, as.numeric) %>% as.data.frame()

window.n <- laply(window.ls, function (win) {
  c(Year=mean(win),
    Close=with(data.use,
       sum(year %in% win & abs(DemMarginGov) < .01 & !is.na(DemMarginGov))))
})
ggplot(as.data.frame(window.n), aes(x=Year, Close)) + geom_point() + geom_smooth()

mypdf("Figure4_GovRDWindows-", width=12, height=4)
(ggplot(rd.win.cast)
 + geom_pointrange(aes(x = as.numeric(as.character(L1)),
                       y = Coef_Conventional, ymin = `CI Lower_Robust`,
                       ymax = `CI Upper_Robust`))
 + geom_hline(yintercept = 0, linetype = 'dotted')
 + ylab('Effect on Policy Change in First Year\n(Robust 95% CI)')
 ## + xlab(paste0("Midpoint of ", (yr.window + 1), "-Year Window"))
 + xlab("Year Range")
 + scale_x_continuous(breaks=seq(1940, 2010, 10),
                      labels=paste0(seq(1940, 2010, 10) - yr.window/2, "-",
                                    seq(1940, 2010, 10) + yr.window/2))
 + theme_bw()
 )
dev.off() 


################################################################################
#### MAIN EFFECT RD ESTIMATES ##################################################
################################################################################

### RD PLOTS
## Equally spaced
par(mfrow=c(2, 2))
ss <- c(1, 2, 5, 10)
optim.bins <- vector('list', 4)
names(optim.bins) <- paste('times', ss)
for (i in seq_along(ss)) {
  s <- ss[i]
  optim.bins[[i]] <-
    rdplot(y=data.use$GovLibD12, x=data.use$DemMarginGov,
           c = 0, binselect='es', scale=s, p=4,
           title=paste0('Evenly Spaced Bins, Scale = ', s),
           x.label='Democratic Share t',
           y.label='Change in Policy Liberalism, t to t+2')
  abline(h=0, lty='dotted')
}

## Bins with equal numbers of units (quantiles)
rdplot(y=data.use$GovLibD12,
       x=data.use$DemMarginGov,
       binselect='qs', scale=1, p=3,
       title='Quantile-Spaced Bins',
       x.label='Democratic Gubernatorial Margin t',
       y.label='Ave. Change in Policy Liberalism, t+1 and t+2',
       y.lim=c(-.11, .11), x.lim=c(-.25, .25))

## Estimates with different "optimal" bandwidths
bw.mserd.tri <- with(data.use,
                     rdbwselect(y=GovLibD12, x=DemMarginGov, bwselect="mserd",
                                kernel='tri'))
(mserd.tri <- with(data.use, rdrobust(y=GovLibD12, x=DemMarginGov,
                                      h=bw.mserd.tri$bws[, "h_l"],
                                      b=bw.mserd.tri$bws[, "b_l"],
                                      all=TRUE, kernel='tri')))
## optimal bandwidth for uniform kernel (simple window)
bw.mserd.uni <- with(data.use,
                     rdbwselect(y=GovLibD12, x=DemMarginGov, bwselect="mserd",
                                kernel='uni'))
(mserd.uni <- with(data.use, rdrobust(y=GovLibD12, x=DemMarginGov,
                                      h=bw.mserd.uni$bws[, "h_l"],
                                      b=bw.mserd.uni$bws[, "b_l"],
                                      all=TRUE, kernel='uni')))
## OLS with clustered SEs in optimal window (linear)
with(subset(data.use, abs(DemMarginGov) < bw.mserd.uni$bws[, "h_l"]),
     robcov(ols(GovLibD12 ~ DemWinGov * DemMarginGov, x=TRUE),
            cluster=abb))
## OLS with clustered SEs in optimal window (quadratic)
## According to what Rocio told me, the confidence interval from this should be
## very close to the robust CI for the linear specification.
with(subset(data.use, abs(DemMarginGov) < bw.mserd.uni$bws[, "h_l"]),
     robcov(ols(GovLibD12 ~ DemWinGov*pol(DemMarginGov, 2), x=TRUE)))
with(subset(data.use, abs(DemMarginGov) < bw.mserd.uni$bws[, "h_l"]),
     robcov(ols(GovLibD12 ~ DemWinGov*pol(DemMarginGov, 2), x=TRUE),
            cluster=abb))


width <- .005
data.use <- mutate(data.use, bin=cut(DemMarginGov, breaks=seq(-1, 1, width)))

bw.mserd.tri1 <- with(data.use,
                       rdbwselect(y=GovLibD1, x=DemMarginGov, bwselect="mserd",
                                  kernel='tri'))
(mserd.tri1 <- with(data.use, rdrobust(y=GovLibD1, x=DemMarginGov,
                                        h=bw.mserd.tri$bws[, "h_l"],
                                        b=bw.mserd.tri$bws[, "b_l"],
                                        all=TRUE, kernel='tri')))

bw.mserd.tri2 <- with(data.use,
                       rdbwselect(y=GovLibD2, x=DemMarginGov, bwselect="mserd",
                                  kernel='tri'))
(mserd.tri2 <- with(data.use, rdrobust(y=GovLibD2, x=DemMarginGov,
                                        h=bw.mserd.tri$bws[, "h_l"],
                                        b=bw.mserd.tri$bws[, "b_l"],
                                        all=TRUE, kernel='tri')))

bw.mserd.tri12 <- with(data.use,
                       rdbwselect(y=GovLibD12, x=DemMarginGov, bwselect="mserd",
                                  kernel='tri'))
(mserd.tri12 <- with(data.use, rdrobust(y=GovLibD12, x=DemMarginGov,
                                        h=bw.mserd.tri$bws[, "h_l"],
                                        b=bw.mserd.tri$bws[, "b_l"],
                                        all=TRUE, kernel='tri')))

bin1.df <- data.frame(bin.mean=tapply(data.use$GovLibD1,
                         data.use$bin, mean, na.rm=TRUE),
                     mid=seq(-1 + width/2, 1 - width/2, width),
                     n=tapply(data.use$GovLibD1, data.use$bin, length))
mypdf("Figure2_GovRD1-", width=8, height=6)
(ggplot(data.use)
 + geom_point(aes(x = DemMarginGov, y = GovLibD1), alpha=.1)
 + geom_smooth(data=subset(data.use, DemMarginGov > 0),
               aes(x = DemMarginGov, y = GovLibD1,
                   weight=tri(DemMarginGov, bw.mserd.tri1$bws[, "h_l"])),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5)
 + geom_smooth(data=subset(data.use, DemMarginGov < 0),
               aes(x = DemMarginGov, y = GovLibD1,
                   weight=tri(DemMarginGov, bw.mserd.tri1$bws[, "h_l"])),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5)
 + geom_point(data=bin1.df, aes(x=mid, y=bin.mean, size=n), pch=1)
 + geom_vline(xintercept=0)
 + geom_hline(yintercept=0, lty='dotted')
 + coord_cartesian(ylim = c(-.15, .15))
 + scale_x_continuous(breaks=seq(-1, 1, .05),
                      limits=c(-bw.mserd.tri1$bws[, "h_l"],
                          bw.mserd.tri1$bws[, "h_l"]),
                      labels=seq(-100, 100, 5))
 + theme_bw()
 + labs(size = '# Obs. in\n0.5% Bin',
        x = "Democratic Margin in Gubernatorial Election (%)",
        y = "Change in Policy Liberalism since Election Year")
 + ggtitle("Governor's First Year in Office")
 + annotate('text', x = 0.001, y = -.025, 0.05, hjust=0, parse=TRUE,
            label = paste('hat(tau)==', round(mserd.tri1$coef['Conventional', ], 3),
                '~(list(', round(mserd.tri1$ci['Robust', 1], 3),
                ',', round(mserd.tri1$ci['Robust', 2], 3), '))'))
 )
dev.off()

